%% Navier-Stokes equation solver
%
%
%

%% Problem settings

clear all
clc
addpath('./source')
global NIt
global NJt
global NIe
global NJe

%% Convergence parameters--------------------------------------------------
convTol = 10^-5;
%--------------------------------------------------------------------------

% Domain dimensions

Lxt  = 12;
Lxe  = 4;       % Longitude before expansion
Lyt  = 2;
Lye  = 1;

% Discretization, pressure finite volumes

% I denotes x direction
% J denotes y direction

elemsPerLenghtUnit = 7;

NIt = Lxt*elemsPerLenghtUnit;
NIe = Lxe*elemsPerLenghtUnit;
NJt = Lyt*elemsPerLenghtUnit;
NJe = Lye*elemsPerLenghtUnit;

% Nodal spacing

h   = Lxt/NIt;

% Constitutive parameters

nu  = 0.01;
rho = 1;

% Boundary conditions


% Time domain and discretization
dt = 0.01;
T0 = 0;
T1 = 100.0;
nSkip = 500;             % skips saving solution
nTimeSteps = (T1-T0)/dt;
nSaveSteps =  (T1-T0)/(dt*nSkip);
nSave      = 0;

dtLim1 = 0.25 * h^2 / nu
dtLim2 = 2*nu/1

if (dt>dtLim1 || dt>dtLim2 )
    fprintf ('stability warning')
end



%% Pre-processing
%  Memory allocation

dimP = NIt*NJt-NIe*(NJt-NJe);
dimU = dimP + NJt;
dimV = dimP + NIt;

% Problem variables
P = zeros(dimP,1);
U = zeros(dimU,1);
V = zeros(dimV,1);

% Auxiliars: midstep variables & Pk for iterative pressure scheme computation
Pk = zeros(dimP,1);
Us = zeros(dimU,1);         % U star
Vs = zeros(dimV,1);


% Saving solution over time
PT = zeros(dimP,nSaveSteps);
UT = zeros(dimU,nSaveSteps);
VT = zeros(dimV,nSaveSteps);

% Auxiliar linear matrix & load vector for pressure implicit equation
M = zeros (dimP,dimP);
l = zeros (dimP,1);


% Dimensions

NIut = NIt+1;
NJut = NJt;
NIue = NIe;
NJue = NJe;

NIvt = NIt;
NJvt = NJt+1;
NIve = NIe;
NJve = NJe+1;

% Mapping spatial points

Xp = PspatialMap( h , Lxt , Lxe , Lye , Lyt , NIt , NJt , NJe , NIe , dimP);
Xu = UspatialMap( h , Lxt , Lxe , Lye , Lyt , NIut , NJut , NJue , NIue , dimU);
Xv = VspatialMap( h , Lxt , Lxe , Lye , Lyt , NIvt , NJvt , NJve , NIve , dimV);

assemblePressureMatrix;

%% Processing - Predictor step
%
%
%
%

% Temporal loop

time = T0;
nstep = 0;
while (time<T1)
    tic
    time = time + dt;
    nstep = nstep+1;
    
    fprintf('Solving timestep %d of %d \n',nstep,nTimeSteps)
    
    
    % Velocity in X
    for i=1:NIut
        for j=1:NJut
            
            if (i<=NIue && j>NJue) continue; end
            
            indexU_cen = UmapOpt (i,j,NIue,NJue,NIut);
            
            U_cen   = U ( indexU_cen );  % center of finite volume
            
            if (j<=NJue)
                
                % Boundary conditions on inlet & outlet
                if (i==1) Us (indexU_cen) = getParabola ( Xp (indexU_cen,2 ) ) ; continue ; end    % Boundary condition inlet
                if (i==NIut) Us (indexU_cen) = Us ( UmapOpt (i-1,j,NIue,NJue,NIut) ) ; continue ; end    % Boundary condition oulet - null normal gradient
                
            else
                
                % Boundary conditions on expansion & outlet
                if (i==NIue+1) Us (indexU_cen) = 0 ; continue ; end    % Boundary condition expansion
                if (i==NIut) Us (indexU_cen) = Us ( UmapOpt (i-1,j,NIue,NJue,NIut) ) ; continue ; end    % Boundary condition oulet - null normal gradient
                
            end
            
            Ax = getAxOpt( i,j,U,V,NIue, NJue, NIut, NJut, NIve, NJve, NIvt, NJvt ,h) ;
            Dx = getDx( i,j,U,V,NJue,NIue ,NJut,h) ;
            
            Us(indexU_cen) = U_cen + dt* (-Ax+nu*Dx);
            
        end
    end
    
    
    % Velocity in Y
    for i=1:NIvt
        for j=1:NJvt
            
            if (i<=NIve && j>NJve) continue; end
            
            indexV_cen = VmapOpt (i  ,j, NIve, NJve, NIvt);
            
            V_cen  = V (indexV_cen );  % center of finite volume
            
            if (j<=NJve)
                % Boundary conditions on top and bot boundaries
                if (j==1) Vs (indexV_cen ) = 0; continue ; end
                if (j==NJve && i<=NIve ) Vs (indexV_cen ) = 0; continue ; end
            else
                % Boundary conditions on top and bot boundaries
                if (j==NJvt ) Vs (indexV_cen ) = 0; continue ; end
            end
            
            Ay  = getAyOpt( i,j,U,V,NIue, NJue, NIut, NJut, NIve, NJve, NIvt, NJvt ,h);
            Dy  = getDy( i,j,U,V,NJve,NIve ,NJvt,h);
            
            Vs(indexV_cen ) = V_cen + dt* (-Ay+nu*Dy);
            
        end
    end
    
    
    %% Processing - pressure step
    %
    
    c2 = rho*h/(dt);
    
    for i=1:NIt
        for j=1:NJt
            
            if (i<=NIe && j>NJe) continue; end
            
            U_right = Us ( UmapOpt (i+1,j,  NIue,NJue,NIut) );
            U_left  = Us ( UmapOpt (i  ,j  ,NIue,NJue,NIut) );
            V_top   = Vs ( VmapOpt (i  ,j+1,NIve, NJve, NIvt) );
            V_bot   = Vs ( VmapOpt (i  ,j  ,NIve, NJve, NIvt) );
            
            if (i~=NIt)
                l ( PmapOpt (i,j,NIe,NJe, NIt) ) = U_right - U_left + V_top - V_bot;
            else
                l ( PmapOpt (i,j,NIe,NJe, NIt) ) = 0 ;
            end
            
        end
    end
    
    l = l*c2;
    
    P = M\l;
    
    
    
    %% Processing - projection step
    
    c1 = (dt/(rho*h));
    
    % Velocity X
    
    for i=1:NIut
        for j=1:NJut
            if (i<=NIue && j>NJue) continue; end
            
            indexU_cen = UmapOpt (i,j,NIue,NJue,NIut);
            
            % Right boundary
            if ( i==NIut ) U ( indexU_cen ) = U ( UmapOpt (i-1,j,NIue,NJue,NIut) ) ; continue ; end
            
            % Left boundary
            if ( i==1 ) U ( indexU_cen ) = Us ( indexU_cen ) ; continue ; end
            if ( i==NIue+1 && j>NJue ) U ( indexU_cen ) = Us ( indexU_cen ) ; continue ; end
            
            
            
            U ( indexU_cen ) = Us ( indexU_cen ) - c1 * ( P ( PmapOpt(i,j,NIe,NJe, NIt) ) - P ( PmapOpt(i-1,j,NIe,NJe, NIt) ) ) ;
        end
    end
    
    % Velocity Y
    
    for i=1:NIvt
        for j=1:NJvt
            if (i<=NIue && j>NJue) continue; end
            
            indexV_cen = VmapOpt(i,j,NIve, NJve, NIvt);
            % Top boundaries
            if ( j==NJvt )
                V(indexV_cen) = Vs(indexV_cen);
                continue ;
            end
            if ( (j==NJve) && (i<=NIve) )
                V(indexV_cen) = Vs(indexV_cen);
                continue ;
            end
            
            % Bot boundary
            if ( j==1 )
                V (indexV_cen) = Vs(indexV_cen);
                continue;
            end
            
            V(indexV_cen) = Vs(indexV_cen) - c1 * (P(PmapOpt(i,j,NIe,NJe, NIt)) - P(PmapOpt(i,j-1,NIe,NJe, NIt)));
        end
    end
    
    % Saving results
    if ( mod(nstep-1,nSkip) == 0 )
        nSave = nSave+1;
        PT (:,nSave)= P;
        UT (:,nSave)= U;
        VT (:,nSave)= V;
        
    end
    
    
    %% convergence check---------------------------------------------------
    if (nstep > 1)
        uDiffSum = 0;
        for i = 1 : dimU;
            uDiff = abs(U(i) - Um1(i));
            if (U(i) ~= 0)
                uDiff = abs( (U(i) - Um1(i)) / U(i) );
            end;
            uDiffSum =  uDiffSum + uDiff;
        end;
        vDiffSum = 0;
        for i = 1 : dimV;
            vDiff = abs(V(i) - Vm1(i));
            if (V(i) ~= 0)
                vDiff = abs( (V(i) - Vm1(i)) ./ V(i) );
            end;
            vDiffSum =  vDiffSum + vDiff;
        end;
        err = uDiffSum/dimU + vDiffSum/dimV;
        errorH(nstep) = err;
        fprintf('Convergence error: %f \n', err);
        if (err <= convTol)
            fprintf('Stationary state reached at iteration %d, convergence error: %f', nstep, err);
            convFlag = true;
            break;
        end;
    end;
    % the previous time step velocities.
    Um1 = U;
    Vm1 = V;
    %----------------------------------------------------------------------
    
    toc
end


%% Post-Process

plotInTime;































